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We develop a nonlinear dynamic theory based on the tensor network state representation that is termed as net- 
work catalyzer dynamics (NCD) to explore the thermodynamic properties of two-dimensional quantum lattice 
models. An imaginary-time-sweep algorithm for implementation of the NCD method is proposed for practical 
numerical simulations. The results obtained by the NCD are disclosed nicely in agreement with the quantum 
Monte Carlo calculations. The quasi-entanglement entropy S , Lyapunov exponent and loop character I p 
are introduced within the dynamic scheme, which are found to display the "nonlocality" near the critical point, 
and can be good characters to detect the thermodynamic phase transitions. 

PACS numbers: 75.10.Jm, 75.40.Mg, 05.30.-d, 02.70.-c 
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The two-dimensional (2D) strongly correlated quantum 
models have triggered broad interest in last decades as they 
always exhibit intriguing and exotic properties such as frac- 
tional quantum Hall effect (H and quantum spin liquid state 
J2l, etc. Accompanied with the boom of quantum information 
science, new theories J^, 0] enable us to describe the criti- 
cal phenomena beyond the traditional paradigm of Landau- 
Ginzburg and renormalization group and access the informa- 
tion of phase transitions by studying the properties of quan- 
tum states without acquiring knowledge of order parameters 
or universality class. Owing to the complexity of many-body 
interactions, reliable and controllable analytical methods for 
these systems are still sparse, and numerical means thus play 
essential roles. Among others, the quantum Monte Carlo 
(QMC), density matrix renormalization group (DMRG) JH1] as 
well as its varian ts tensor network state (TNS) based 

algorithms JI- IT, 1511 and so on, have achieved great success 
1 12l Il3ll . However, QMC is not applicable to the frustrated 
spin systems and Hubbard model away from the half-filling 
because of the "negative sign" problem, and DMRG is re- 
markably accurate and efficient in one dimension but has great 
costs for 2D systems of large size. Therefore, to develop new 
theories and efficient algorithms for correlated quantum lattice 
systems is highly encouraged. 

In this Letter, we develop a nonlinear dynamic theory with 
the TNS representation of the density operator, termed as the 
network catalyzer dynamics (NCD), for exploring the ther- 
modynamic properties of the infinite 2D quantum spin lattice 
models. The partition function as well as the averages of ob- 
servable operators can be expressed in a simple form with the 
cell tensor and the corresponding catalyzers that can be ob- 
tained through a dynamic way. Within the framework of the 
NCD, the quasi-entanglement entropy S defined by the trans- 
fer matrix of the partition function, Lyapunov exponent I lya 
that quantifies the convergent properties and loop character 
ji°°p f or the loop dependence are introduced to describe the 
properties of the TNS. An imaginary-time-sweep algorithm 
that is free from the negative sign problem for implementation 
of the NCD method is proposed. The results calculated by the 
NCD for the spin- 1/2 Heisenberg antiferromagnet (HAF) on 



honeycomb lattice are nicely compared with the QMC simu- 
lations, showing the efficiency and accuracy of such a method. 
The quantities S, l' ya and l' oop are revealed to be able to de- 
tect possible thermodynamic phase transitions of the system, 
which are all characterized by sharp peaks at critical points, 
as manifested by the spin- 1/2 anisotropic honeycomb HAF 
model as an example. 

Tensor network state representation for the density 
operator. — Suppose that the Hamiltonian can be written as 
H = H'i, where is the local Hamiltonian of two con- 
nected spins at zth and y'th lattice sites. Without losing gener- 
ality, we take the infinite honeycomb lattice as the prototype 
here. Let us begin with introducing the local evolution oper- 
ator U ij = e- jH " = Zptpjp'p' Up^PjPpiPjXPiP'jl where T is 
the infinitesimal imaginary time slice, \pi) is the local basis 
of the zth spin and p, is called the physical bond. By using 
the Trotter-Suzuki decomposition IU4II . we can write the den- 
sity operator at inverse temperature j3 as p(J3) = [Y\{ij) U'i] K 
with f3 = Kt. Define the evolution tensors G L and G R as 

X, where C/^ p , g 



G , = U L , 

PiP,g PiP r g 



'*< and G Ug = U U, 



(U R p , g ) contains the left (right) singular vectors and A g con- 
tains the singular values of the matrix Uip lP >)ipjp>.). We call 
g the geometrical bond. Then the partition function is rep- 
resented as the contraction of a closed tensor network (TN), 
which is a three-dimensional brick-wall structure formed by 
U L 's and U R 's. Obviously, it is impossible to complete such 
a contraction exactly. One way to do this is first to con- 
struct the tensor product density operator (TPDO) H5Y at in- 
verse temperature r as p(r) = gTr{\[ i=odd A' Ylj =even W) with 
A'' , = y D ",y-G L „ G L G L „, , , B j , 

PiPpgagbgc '-'Pi Pi PiP",ga P"P, ,gb P'"Pj,gc' PjPpgagbgc 

y,y<,y« G R „ G R G R ,„ , , and then to evolve the TPDO 

^PjPj Pjp'j,ga P"p"',gb P"Pj,gc' 

as A' „. = Yn'A' , G L , „ , with the composite in- 

PiP,,g,igbgc ^Pi PiP\,gagbgc PjP",g„ ^ 

dex g a = {ga,g' a ) while keeping its form unchanged (the con- 
tractions of B-' and G R are similar). The gTr means sum- 
ming over all shared geometrical bonds. During the evolution, 
the geometrical bonds will be enlarged exponentially with the 
contractions, and proper approximations are needed to bound 
their dimensions. 



Network catalyzer dynamics. — Denote f3 as the inverse 
temperature of the present TPDO p(J3), and (J3 + /?') is the 
targeted inverse temperature. The partition function may be 
rewritten as Z(J3 + j3') = Tr[p(J3)p(J3')]. Suppose that we also 
have the TPDO representation of p(J3') that is comprised of 
A's and B's. Consequently, 



Z(p+p') = P Tr[gTr(Y] A' ] ] B VX [~[ A? 

i—odd j-even 



i-odd 



(i) 



j=even 



where pTr stands for summing over all shared physical bonds, 



and <^ 



are the Dirichlet matrices with elements S p >p t and 



6p' kPk , respectively. 

Now we define the catalyzable closed TN of the partition 
function with the following two restrictions: (a) The TN is 
formed only by one inequivalent tensor that is called the cell 
tensor and denoted as T"'". The cell tensor may be chosen 
as a small cluster that consists of several properly selected 
original tensors. Straightforwardly, a new larger cell tensor 
can be constructed by a cluster of cell tensors. We denote 
T cc " (1) as the smallest cell tensor and T re " (y) as the cell tensor 
that is made up of y T ee "W's, where y is called the size of T' e ". 
(b) If the t/th bond of a T ce " is connected with the vth bond 
of its adjacent T ce " with u + v, then under the permutation 
operation between the t/th and vth bonds, the T ceU is required 
invariant, say T ce ' 1 , = T ce f . . We call the bonds f„ and 
£, v an index pair. When u = v, g u forms a pair with itself. 

For the spin- 1/2 Heisenberg model on the infinite honey- 
comb lattice, we can rewrite Eq. ([TJ into the catalyzable form 



as Z = gTr(Ui T '') with T 



n cell ,i 



3l c £ , and 

4-4 <5 



3\ 



= V A' , Al dK-Sl , (2) 

flftfi £_i PiP r glg2g3 PiP r glg2g3 PjPi P,P:' V ' 



- E 



B } , Si 6 J ,. Si , (3) 

PiPj,g4g5g6 PjPj,gig5g6 PjPj PjPj' V ' 
PjP'jPjPj 

where the composite bond = (g a ,g a ) [Fig. [U(a)]. The par- 
tition function is now represented as a square closed TN [Fig. 
Q](b)]. In T cc ", £i (£2) forms an index pair with £4 (£g). Owing 
to the symmetry of the model, the permutation invariance is 
satisfied. 

Now we introduce the catalyzers x a (a = l,2,...,J with 
/ the number of index pairs) that are a set of unit-norm col- 
umn vectors defined on a catalyzable TN formed by the nth- 
order T"'" with following conditions: (a) Recursive condition. 
The catalyzers x ff are the solutions of the set of self-consistent 



fcell Y ~ 



equations rx." = 2f 

with T a positive constant for any a = 1, . . . , n. Henceforth, 
a a = ctb if ath and feth bonds form a pair. The T ce " is consid- 
ered as a nonlinear mapping T ce " : Y\®a*b^ a " ~> V ff " with 
V a ° the vector space of x" n . To describe the convergence of 
the mapping the Lyapunov exponent I lya can be defined as 



fy = lim VllnV l —— _, V -}/0, (4) 

8=1 a=\ 




FIG. 1: (Color online) (a) The closed TN that is comprised of two 
inequivalent tensors and S, which represents the partition func- 
tion of the spin- 1/2 Heisenberg antiferromagnet on the honeycomb 
lattice, (b) The catalyzable TN formed by T™" is obtained by con- 
tracting a common bond between 3\ and S. (c) One possible config- 
uration of the defective TN in which certain T""'s are replaced by 
the rank-1 tensor T"'" to forbidden the loops formed only by T CE "'s. 



with e a ° the infinitesimal random vector and T e = T[T(. . .)]. 
(b) Optimization condition. The least-squares cost function 

D = JZ^JT™» 2 ... -n fl ^;] 2 is minimized 03]. An 
equivalent representation of the optimization condition is that 
the inner product function f = 2fif 2 ... T^l 1 x^x^ • ■ ■ > 
is maximized. From the optimization condition, the x ff 's 
form the rank-1 decomposition of the T eefl , namely T ceH ^ 
Tlig)a xC ""- m what follows, we shall show that the infinite 
contraction for the partition function as well as observables 
can be dramatically simplified in the framework mentioned 
above by approximating the catalyzable TN with defects. The 
loop character I loop is thus proposed to display the validity of 
such an approximation. 

The defective TN is obtained by replacing minimal number 
of T Cf, "'s in the catalyzable TN with the rank-1 tensor T ce " = 
il(g)fl xl? "> so mat there exist no loops in the TN formed only 
by T ce " [Fig. Q](c)]. The f ce " is called a defect. Under this 
circumstance, the partition function Z of the prototype model 
in Eq. (HJ is simplified as 



n\e a °\ 



Z = Y N - l gTr(T e "t ce ") = t N - l gTr(T e " \\ = f", (5) 



with N — > 00 the number of the T celh s. Equation © can be 
proved by reversing the infinite contraction procedure into an 
infinite growing one. From Eq. (0, the defective TN can be 
naturally rebuilt with minimal number of T cp "'s. The recur- 
sive and optimization conditions are employed repeatedly to 
grow the TN and construct the T ee " with catalyzers, respec- 
tively. We illustrate a part of the growing procedure shown in 
Fig. |2](a). It should be cautious that the defect reconstructed 
is Tfjj* c , and the permutation invariance (restriction (b) of 
the catalyzable TN) is needed to ensure this defect is the rank- 
1 tensor of the original Ti el i , , . One may note that although 
there are many possible configurations for which cell tensors 
will be replaced by the defects, Eq. (0 holds so long as the 
loops of pure T^'s are forbidden. Alternatively, one can be- 
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FIG. 2: (Color online) (a) A part of the growing procedure of the 
defective TN beginning from Eq. {5). The recursive condition is used 
to grow the TN and the optimization condition is used to construct 
the rank-1 tensor T re " with x"''s. (b) One possible construction of 
the partition function with loops of X ceH(1 ''s, 



gin with the catalyzers and a cluster of the cell tensors that 
already contains the loops of T' e "'s [Fig. 12(b)], and grow the 
rest, but in this case, Eq. (0 should be modified accordingly. 

The thermal average of an operator 0' J '"•* can be similarly 
simplified within the NCD scheme. It can be noticed that, 
except the tensors that share the physical space with O lJ the 
rest of the TN are exactly the same as those of the partition 
function. We formally denote the cluster formed by tensors 
sharing the physical space with O'j--* as Q>p i ...p k p>....p' tgl ... gN . The 
average 0' J ~ k ) can be obtained by 

N 

(d ij - k )=pTr[d ij - k gTr(<S> ] ] x a ')]/Z 

N 

= I Kin-.,: X<<V.,> | |vM//. 

Pt-Pk gu-gN a=\ 

P'i-P'k 

(6) 

The approximation can be improved by increasing the size 
of the cell tensor. For the present HAF model, we can con- 
struct T cel1 through T 1 " and T- 1 -: 

j,cell(y) _ yTn jl7i jl72 

(a l a])(b 3 b i )(c 2 c 4 )(did 2 ) /—I a\biC\d\ a 1 b 1 c 1 d 1 a 3 b 3 adi 

a2a^b\b2C\c^d^di r 

T S i o i d i S ^ S bid3 GavcSbvh > (7) 

with y = 2{y\ + yi). T T and T- 1 are initiated as T ce " (1) 
and increased as T 1 ^,^ = T%cdKt% 6 bd' and 
T WMcc')d> = ^Vd^atd T tv% 6 Vd, respectively. Note that 
this is not the only way to increase T ce "'s size, but is the way 
that can preserve the symmetry of T ce " with truncations at a 
reasonable cost. To see how the increase of y improves the 
approximation, we introduce the "transfer matrix" of the par- 
tition function, M(y) = g T r[T cen{y) x a - y x a - y }. Thus we can 
simplify Z at large y as Z lim;^ M f 'TrlAiiy)'], where the 
catalyzer x ab < d > ,r is the dominant eigenvector of M(y). On the 
other hand, as x a '' u "' 7 can be written in the form of a matrix 

jr— i 

product state (MPS) 117!] with {y\ + yi) "physical" bonds, we 



may define the quasi-entanglement entropy (QEE) of the MPS 

as 

s^-Yj^T^fm^fi (8) 

with fi^' 7 the normalized singular spectrum of the matrix 

m„y , a„,y, 
x b i b i y^didf- 

In addition, one can see that the defective TN with 
T «>«(r+i)' s and f««(y+i)» s conta i ns more loops of T ceUm than 

that with T cefl W's and f cc "W's (recall that in this defective TN, 
loops of pure T cc " (r) 's, not T^^'s, are forbidden). In this 
way, we can choose the catalyzers of the unchanged bonds b 
and d to define the loop character I loop as 

jloop(y) _ 2 _ ( x ai M) ,y+l^T x a m ,y ^ 

with x a " 7 the catalyzers of T ce " (r) . Such a character I locp 
quantifies the loop dependence of the TPDO. It is worth men- 
tioning that while I lya delineates the properties of a certain 
mapping, I loop describes the difference between two mappings 
with different cell tensor sizes. Specifically, when r ™ con- 
verges to zero as y increases, the effects from larger loops are 
negligible. So, I Io "pW can be used to monitor the error brought 
by the defective TN approximation. 

The dimensions of the cell tensor's bonds a and c are en- 
larged during the growing procedure while those of the bonds 
b and d are unchanged. The truncation to bound the dimension 
is the same as the way to bound the dimension of the geomet- 
rical bond of the TPDO during the imaginary time evolution, 
as will be introduced later. 

The power al gori thm for catalyzers. — We employ the 
power algorithm [16] to calculate the catalyzers of T cc ". The 
x' ,; is renewed alternatingly until its convergence by 

rj-cell^mii+Y) x «,_i(r+l) x » i+ i(r) j > ^anit+l) ^q-j 

with t the iteration times. The Lyapunov exponent I lya can be 
utilized to describe the property of convergence for a given 
T cc ". In particular, for a given convergence with smaller I lya 
less iteration times are needed, and vice versa. Notice that 
one iteration step for the catalyzers with the mapping ( fTOt is 
actually equivalent to growing one T ee " in the defective TN 
picture. In other words, the longer time to converge is equiv- 
alent to the contraction of a larger defective TN. Thus, an in- 
teresting meaning of I lya in the NCD is that it indicates the 
"spatial-dependence" of the wave function with its local ten- 
sor representation. In our calculations, I lya keeps negative for 
all parameter ranges. For I lya being zero or positive, the sys- 
tem might have possible bifurcations or chaos. 

The imaginary-time-sweep algorithm. — For the TPDO at 
finite temperature, we propose the imaginary-time-sweep al- 
gorithm (ITSA) based on the NCD scheme. First, we intro- 
duce the truncation principle for the enlarged bonds. Denote 
the bond to be truncated as g with its dimension x an d define 
the matrix M such that Z = £ jr hi'^W = ^(M), which 
is obtained by simply contracting all shared bonds in Eq.(Q]i 
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except g. According to the linear algebra, the best truncation 
of the bond g can be reached using the singular value decom- 
position (S VD) on M, say M PAQ r , where only x largest 
singular values and the corresponding left and right singular 
vectors are kept. Here x is the preset dimension cut-off. The 
truncation error can be controlled by 



a= X +\ 



(ID 



Redefining the matrix as M = VAP r Q VX, we have Z = 
7>(M), and then, use the SVD again to decompose M as M = 
PAQ r . The best^ xx truncation matrices for bonds g and g' 
are obtained by 



P = PA- 1/2 P VX, Q = QA 1/2 Q VX, 



(12) 



such that 6 = PQ T (6 is the Dirichlet matrix). Thus, as long 
as we have M, we can get the best truncation matrices. Obvi- 
ously, the difficulty in obtaining M is as big as in calculating 
Z itself. To simplify the calculation of M, we may adopt the 
defective TN scheme of the NCD in the way of 



M„ 



f N 



x\ ,X J , „ 

Aatedfg'l ^Mgjg'J S iffg'jg'J 4„ X l X i /f ' 



(13) 



where we write selectively the composite index into (ghg'j) 
for clarity. The first line of the r.h.s. of Eq. $13[ holds when 
the enlarged bond is one of T ce "'s bonds, and the second line 
holds when the enlarged bond is contracted in the construction 
of T cc ". Equation ( [TBI can be proved by infinitely growing the 
defective TN, similar to that for calculating the observables. 

There are two main parts of the ITSA (Fig. |3), initializa- 
tion and sweeping procedure. In the initialization part, each 
TPDO at the inverse temperature /? < /? is initialized is the 
targeted temperature). The M is calculated with the partition 
function Z(2/3) constructed by the TPDO p(J3) and its copy, 
where /3 that begins from r is the inverse temperature of the 
present TPDO. In the sweeping part, the M for each trunca- 
tion is obtained from the partition function Z(J3), i.e. when 
the TPDO is evolved to denoted as p(J3), Zifi) can be con- 
structed with p(J3) and the pre-obtained p0 - p). If it is the 
first sweep for this targeted temperature fl, we can construct 
Z(j6) using p(j3-/3) saved in the initialization, otherwise using 
p(J3 - /}) obtained in the last sweeping step. 

The errors in the ITSA come from three parts: the error 
from the Trotter- Suzuki decomposition, the truncation error 
and the error caused by introducing the defects (rank-1 ten- 
sors). The first part is controlled by r, and the second is 
controlled by s [Eq. (fTTTil. The third part affects the cal- 
culations of transformation matrices and observables, and a 
possible choice to control it is the cost function D (that is in 
fact the error of the single tensor rank-1 approximation). But 
considering that the contribution from the "residual" tensor 
T ce " - T ce " may vanish along with the infinite tensor product, 



| Step 1: Obtain G L and G R and calculate the TPDO p(P=t).| 



Step 2: Evolve the TPDO p(P) with G L and G R to p'(P). Use p'(P) and 
its copy to construct Z(2P) and calculate the transformation matrices. 




Truncate the enlarged bonds of p' and save the TPDO. a. 



( 1 ) Stop when the TPDO's at all targeted 
temperatures are obtained, or (2) go to step 2. 



Step 4: Load G L and G R and the initial TPDO p(P=x) . 4 




Step 5: Evolve p(P) with G L and G R to p'(P), use p'(P) 
and the pre-saved p(P-P) to construct the partition 
function Z(p) and calculate the transformation matrices. 



Step 6: Truncate the enlarged bonds of p'(P) and save the TPDO. 



(1) Stop when the present temperature reaches p or (2) go back to step 5. 



Stop when the TPDO p(P) convergences or go back to step 4. 



FIG. 3: (Color online) The sketch of the ITSA algorithm. 



the loop character I loop is a more solid choice to control the 
defect-caused error. 

Applications to the spin-1/2 Heisenberg antiferromagnet on 
honeycomb lattice. — To show the efficiency and accuracy of 
the NCD approach, we take the spin-1/2 Heisenberg antifer- 
romagnet on honeycomb lattice with nearest neighbor inter- 
actions as an example, whose local Hamiltonian reads W ! = 
A[§Kx)§Kx) + §i(y)§M] +S i{z) S jiz \ where A characterizes the 
anisotropy of exchange interactions. First, we calculate the 
energy per site E = Tr(Hjjp) at A = 0.5 and 1 to testify the va- 
lidity of NCD [Fig. H(a)]. The translation invariance has been 
used. The results are compared with QMC simulations, in 
which different dimension cut-offs^ are used in the NCD cal- 
culations. We use the tensor cluster in Fig. |2](b) to calculate 
the truncation matrices, where the truncation error e is found 
around 10 5 . The iteration in the power algorithm is stopped 
when the difference of the catalyzers at t and t + 1 steps [Eq. 
([TOi ll is less than 10~ 14 . In calculating the observables within 
the NCD, the loop character l'" op is required under 10~ 7 and 
the cell tensor size is increased until the difference of between 
the observables with size y and y + 2 is less than 10~ 6 . We set 
the lattice size as 64 x 64 for QMC calculations and keep the 
errors around 1CT 5 . For the characters of the TNS, we fix the 
cell tensor size y = 500 and set = 300 [Eq. ©] for l' ya . 

It is known that in the present system a thermodynamic 
phase transition (TPT) may occur in the presence of the spin 
anisotropy A, and the critical temperature can be determined 
by the divergent peak of specific heat C = -/} 2 dEld[3. Our 
calculation of C shows that the TPT occurs at T c = 0.345 
with A = 0.5 [Fig. gl(b)] for both NCD and QMC calcula- 
tions. We also calculated the QEE [Eq. ©], where the similar 
behavior near the critical point is observed [Fig. 0(a)]. The 
sharp peak of S appears at T s = 0.3610, close to the value 
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FIG. 4: (Color online) (a) Energy per site E obtained by NCD and 
QMC at A = 0.5 and 1. (b) The staggered magnetization per site 
m z and the specific heat C as functions of temperature T for A = 0.5. 
The specific heat indicates a thermodynamic phase transition at T c = 
0.345. 
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FIG. 5: (Color online) (a) The temperature dependence of the QEE S, 
in which the critical temperature T s = 0.3610. (b) The temperature 
dependence of the Lyapunov exponent and the loop character, which 
show sharp peaks at T lya = 0.3636 and T 1oo p = 0.3597, respectively. 



obtained from the specific heat. The QEE vanishes (about 
10~ 3 ~ 10~ 4 ) when the temperature is away from the critical 
vicinity. 

We studied the ^-dependence of the Lyapunov exponent 
I'y [Eq. ©] and the loop character/' 00 '' [Eq. ©]. The re- 
sults [Fig. |5](b)] show that the characters reach maximum at 
T l y = 0.3636 and r /o °'' = 0.3597, again close to the critical 
temperatures obtained from the specific heat and the QEE. I ly " 
shows that the TPDO near the critical point has the property 
of considerably strong spatial-dependence, suggesting that the 
TPDO could be well approximated with a defective TN of 
great size. From /' oop one may observe that the TPDO bears 
high loop dependence and the effects from the loops of cell 
tensors can be ignored only when the size of cell tensor y is 
large enough. Meanwhile, noticing that the QEE versus T 
behaves differently in different phases, showing promisingly 
a prospective development of NCD, that is, the phases can 
be classified by the behavior of the QEE. When the temper- 
ature is away from the critical point, the whole TN can be 
approximated by a relatively small defective TN, where the 
effects from the loops of cell tensors can be ignored at rela- 
tively small size. It should be remarked that as it is difficult 



and I loop will be coincident. 

In conclusion, we developed the NCD for exploring the 
thermodynamic properties of 2D quantum lattice models, and 
proposed the ITSA that is free from the negative sign problem 
for the numerical realization of NCD. The calculated results 
for the spin- 1/2 HAF on honeycomb lattice as an example are 
found consistent well with the QMC simulations. New char- 
acters such as the quasi-entanglement entropy, Lyapunov ex- 
ponent and loop character are suggested to describe properties 
of the thermal states and detect possible thermodynamic phase 
transitions of the system. The possible extension of the NCD 
to other quantum lattice systems can be done in a straightfor- 
ward way. 

The authors are indebted to W. Li, X. Yan, Y. Zhao, 
Z. C. Wang and Q. R. Zheng for stimulating discussions. 
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for a TNS to approximate the state in the critical vicinity II 1 811 . 
it is expected that with sufficient large the critical tempera- 
tures determined by the specific heat and the characters S , I lya 
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